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ABSTRACT 



1 INTRODUCTION 



We discuss results of global three-dimensional (3D) magnetohydrodynamic (MHD) 
simulations of accretion on to a rotating magnetized star with a tilted dipole magnetic 
field, where the accretion is driven by the magneto-rotational instability (MRI). The 
simulations show that MRI-driven turbulence develops in the disc, and angular mo- 
mentum is transported outwards due primarily to the magnetic stress. The turbulent 
flow is strongly inhomogeneous and the densest matter is in azimuthally-stretched tur- 
bulent cells. We investigate two regimes of accretion: a magnetospheric regime and a 
boundary layer (BL) regime. In the magnetospheric regime, the magnetic field of the 
star is dynamically important: the accretion disc is truncated by the star's magnetic 
field within a few stellar radii from the star's surface, and matter flows to the star in 
funnel streams. The funnel streams flow towards the south and north magnetic poles 
but are not equal due to the inhomogeneity of the flow The hot spots on the stellar 
surface are not sjTnmetric as well. In the BL regime, the magnetic field of the star 
is dynamically unimportant, and matter accretes to the surface of the star through 
the boundary layer The magnetic field in the inner disc is strongly amplified by the 
shear of the accretion flow, and the matter and magnetic stresses become comparable. 
Accreting matter forms a belt-shaped hot region on the surface of the star The belt 
has inhomogeneous density distribution which varies in time due to variable accretion 
rate. The peaks in the variability curve are associated with accretion of individual tur- 
bulent cells. They show 20-50% density amplifications at periods ^5 — 10 dynamical 
time-scales at the surface of the star Spiral waves in the disc are excited in both, mag- 
netospheric and boundary layer regimes of accretion. Results of simulations can be 
applied to classical T Tauri stars, accreting brown dwarfs, millisecond pulsars, dwarf 
novae cataclysmic variables, and other stars with magnetospheres smaller than several 
stellar radii. 



boundary layer regime (e.g., Popham & Narayan||1995 Inog 



Different types of stars have dynamically important magnetic 
fields. These include young, classical T Tauri stars (here- 
after - CTTS; e.g., [Bouvier, et a l. 2007), young accreting 
brown dwa rfs (e.g.,|Mohanty et al.,2005j, ac creting neutron 



stars (e.g., |Lewin et al.||1995| [Van der Klis||2000| l, and dif- 



ferent t5^es of cataclysmic variables, including dwarf novae 
and intermediate polars (e.g., Hellier 20011 | Warner|[2003[ l. 
In such stars, the accretion disc is truncated at the disc- 
magnetosphere boundary, and the magnetic field governs the 



matter flow (e.g.,|Pringle & Rees|1972t|Ghosh & Lamb 


1979 


K6nigri991 
Campbell 2( 


'Collier Cameron & Campbell' 1993 


Wang 


1995 


nOl. We refer to this regime as t 


le magneto- 



am ov & SunyaevI 19991 |Piro & Bildsten|2004||Fisker & Balsara 
(2000^ 

The physics of the disc-magnetosphere interaction de- 
pends on both the properties of the magnetized star and the 
properties of the accretion disc. Accretion in the disc can be 
smooth as modelled by including an a— viscosity ( [Shakura~&| 
|Sunyaev|1973P in the MHD equations. Or, it can be strongly 
turbulent owing to the growth of the magneto-rotational in- 
stability (Velik hovl 19591 |Chandrasekhar| 1960} [Balbus & Haw-| 
|ky[r991,,1998) . 

Interaction of a magnetized star with a— type discs has 
been investigated in a number of axisymmetric simulations 



spheric regime of accretion. On the other hand, many accret- 
ing neutron stars, white dwarfs, and other stars are expected 
to have relatively weak, dynamically-unimportant magnetic 
field, and accrete in the equatorial region of the star in the 



(e.g., ,Miller & Stone 1997 


1 Romanova et al. 


2002( |Long| 


|et al.| 2005 |Bessolaz et al. 


20081, and in global 3D simu- 


lations (e.g., [Romanova et al."2003 


2004 2008 Kulkarnil 


|& Romanovaf|2005 , 2008; Long et a 


. 2007 2008 20111. 



These simulations confirmed many properties of the disc- 
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magnetosphere interaction predicted by theory. They also re- 
vealed new features, connected, for example, with the non- 
stationarity of the disc-magnetosphere interaction and the 
possibility of outflows (e.g., Lovelace et al. 1995[ Goodson 
et al.l|1997] |1999t [ Romano va et al.||2005[ [Ustyugova et al. 
2006 ( 'Roma nova et al.|2009 l, and possible important role of 
the interchange instability l |Romanova et a l. 2008 ; Ku lkarni &| 
Romanova|20 08 2009 ). Interaction of a magnetized star with 
a turbulent accretion disc has not been studied as much. MRI- 
driven turbulence has been extensively studied in axisymmet- 
ric and local/ global 3D MHD simulations (e.g.,'Hawley et al.| 
1995; Brandenburg et al. 1995; Stone et al. 1996; Armi tage , 
1998; Ha wley et al . 2001 ; Stone & Pringle 2001 ; Bec kwith7et| 
al.„2009||201 It [Flock, et al..2011,, Simon et al, 20111 ). How- 



ever, in all these simulations, the central object is a black hole. 

Recently, we performed the first axisymmetric (2.5D) 
MHD simulations of the MRI-driven accretion on to a mag- 
netized star ( Romanova et al.]|2011ai ). We observed differ- 
ent phenomena connected with turbulent nature of the disc. 
However, global 3D simulations are required for more realistic 
description of MRI-driven turbulence and magnetized stars. 

In this paper, we present results from global 3D MHD 
simulations of MRI-driven, turbulent accretion on to a magne- 
tized star with a tilted dipole magnetic field. It is necessary to 
perform the simulations in the full three-dimensional space, 
including the full range of the azimuthal angle (0 — 2it) . Fur- 
thermore, we need to have sufficiently high grid resolution in 
the disc for the magnetic turbulence to develop and be sus- 
tained, and we need to resolve the magnetosphere of the star, 
which has strong gradients of the dipole magnetic field. Be- 
cause of these requirements, significant computing time is re- 
quired, and we performed simulations only during a restricted 
period of time. 

First, we show results of MRI-driven accretion on to stars 
with dynamically unimportant magnetic field, accreting in the 
BL regime. In this case, we were able to obtain relatively long 
runs and investigated MRI-driven accretion in the disc. We 
also studied the properties of the inner regions of the disc, the 
disc-star boundary, and variability associated with the turbu- 
lence. Next, we show results for the magnetospheric accretion 
where the disc matter is truncated by the magnetosphere. In 
this case, we concentrated on the disc-magnetosphere inter- 
action and accretion through the funnel flow. 

In §2 we describe the problem setup. In §3 and §4 we 
show results of simulations for accretion in the boundary 
layer and magnetospheric regimes, respectively. In §5 we dis- 
cuss the different stresses in the disc. §6 gives the conclusions 
from this work. 



2 PROBLEM SETUP 

2.1 Magnetosphere and magnetospheric radius 

If the magnetic field of the star is relatively strong (in the 
sense of being dynamically important]]^ then the accretion 
disc is truncated by the stellar magnetosphere at some ra- 
dius r,n, which is called the magnetospheric radius. Different 
physical arguments were proposed to define this radius (e.g.. 



^ Here, and below in the paper, we use terms 'strong' and 'weak' for 
dynamically important and dynamically unimportant fields. 




Figure 1. The "cubed sphere" grid used in our simulations. The grid 
is compressed towards the equatorial plane. A low resolution case is 
shown here to make the grid visible. 



[Pringle & Rees|1972[|Ghosh & Lamb|1979p . For example, this 
radius can be derived from the balance between the magnetic 
pressure of the mag netosphere and the ram pressure of the ac 



creting matter (e.g.. Lamb et al.|l973 1: B^/Stt ^ pv'' (where 



B is the local value of the magnetic field, p and v are den- 
sity and total velocity in the disc) . Then, the magnetospheric 
radius can be derived in the form, which is often used in the- 
oretical analysis (e.g, [Lamb et al.|1973][K6nigl|19"9T] l: 

(1) 



-2/7 



GM 



-1/7 



where p — Bi,Rl is the magnetic moment of the star, M is the 
accretion rate from the disc, and coefficient fc ^ 0.5 l [Long et[ 
,al.,2005) . 

In simulations, however, we do not know the accretion 
rate in advance, and it varies with time. Hence, it is better to 
find the magnetospheric radius directly from the balance of 
stresses. We introduce 



Pi = 



P 



B2 



(2) 



where, p is the gas pressure, is the azimuthal component 
of the magnetic field. Here, /3 is the commonly used plasma 
parameter. The modified plasma parameter /3i takes into ac- 
count the matter stress associated the plasma flow in the disc 
and it is often more relevant in astrophysical situations than 
/3 parameter (Romanova et al. 2002). In accretion discs, the 
magnitude flow velocity is approximately the azimuthal ve- 
locity, V ^ Vtf,, and hence the matter stress has also a meaning 
of the ram pressure, and the condition /3i = 1 can be also in- 
terpreted as a balance of total pressure at the disc-star bound- 
ary. We observed from simulations that this condition always 
gives the boundary between the disc and the magnetosphere. 
We use this condition to find r-m from simulations. Typically 
we have /? ^ 1 at radii larger than rm, which often shows the 
place where the funnel streams start flowing from the disc to 
the star (see also Bessolaz et al. 2008; Campbell 2010) . p| 

Modeling of accreting magnetized stars is technically 
challenging, because in regions of the strongest magnetic 
field, the matter density is low, and this requires a small time- 
step in simulations. In a— type discs, we were able to model 
relatively large magnetospheres with r,n ^ WRi, (e.g., [Ro-[ 
[manova et al.|2006| ). However, in the current simulations of 



Bessolaz et al.[ l [2008 t compared several approaches for derivation 



of r^n and showed that this radius does not differ much. This is be- 
cause the magnetic pressure ~ ~ r"^ is a steep function of r. 
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Figure 2. The initial distribution of density (color background) and magnetic field lines in the xz— slice. The left-hand panel shows the case 
of a dynamically unimportant and slightly tilted stellar magnetic field (Bl = 0.001, = 5°). The two right-hand panels show the cases of a 
dynamically important stellar magnetic field (B'^ = 10) at high (6 = 30°) and low (0 = 2°) tilts of the dipole field. The magnetic field in the 
0.01, respectively. 
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Figure 3. Slices of the density distribution (color background) and sample field lines in the case of a star with dynamically unimportant magnetic 
field (B', = 0.001, = 5°). The left-hand panels show the xz and yz slices, while the right-hand panel shows an equatorial, xy, slice. 



turbulent discs, we model somewhat smaller magnetospheres 
with r„i ^ 37?*, because the duration of simulations strongly 
increases with the size of the magnetosphere. There are sev- 
eral types of stars where such magnetospheres are expected 
(e.g., CTTSs, accreting millisecond pulsars, dwarf novae, ac- 
creting brown dwarfs, and other stars) . 

In case of much weaker, dynamically unimportant field, 
the star's magnetosphere does not truncate the disc. That is, 
the magnetospheric radius is equal or smaller than the ra- 
dius of the star, < i?,. We use this case for investiga- 
tion of the MRI-driven accretion in the disc and for analy- 
sis of processes at the disc-star boundary in the presence of 
the magnetic field. These simulations are faster, and they also 
correspond to an important astrophysical regime of accretion 
through the boundary layer. 

2.2 Rotation of the star and corotation radius 

Another important parameter of the problem is the corotation 
radius, rcor = (GA/*/f2*)^/^, where Keplerian angular veloc- 
ity in the disc equals the angular velocity of the star, . The 
result of the disc-magnetosphere interaction depends on the 



ratio between r,n and rcor- In case of slow rotation, rcor fm, 
the magnetosphere spins the inner parts of the disc down, the 
gravitational force is larger than the centrifugal force, and ac- 
cretion is favorable. In the opposite case of rapidly rotating 
star, rcor < r™,, the magnetosphere transfers its angular mo- 
mentum to the matter of the inner disc, and a star is in the 
'propeller' regime, where accretion is suppressed (e.g.,|Illari- 



onov & Sunyaev|1975[ Lovelace et al.|1999}|Romanova et al. 



2005 Ustjoigova et al. 2006 1. In this paper we consider only 
the case of slowly rotating stars, which is favorable for accre- 
tion. We chose rather large corotation radius, rcor = 10-R*, to 
be sure that we are always in the accretion regime from the 
beginning of the simulations]^ The ratio rcor/^^m ~ 3 is larger 
than that expected in the rotational equilibrium stat^ where 
rcor ~ (1.2 — 1.6)rm (the coefficient in front of r^ depends 
on parameters of the model l [Long et al.|2005p ). However, our 

^ Initially, the inner disc radius is located at r = lOi?*, and this fact 
determines our choice for rcor- Smaller values of rcor would be also 
sufficient. 

* In the rotational equilibrium state, the torques on the star associ- 
ated with spinning-up and spinning-down are balanced. 
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earlier simulations of a— discs, performed at different ratios 
of Tm/rcor, show that the result is very similar in cases where 
Tcor is larger, or much larger than Vm- In both cases, the gravi- 
tational force dominates and the accretion through the funnel 
flow has been observed (e.g., |Romanova et al.|2002(|Long et| 
|al.|2005| |. 

2.3 MRI-driven turbulence 

The turbulence in the disc is initiated and supported by 
the magneto-rotational instability. Here, we briefly summa- 
rize the condition for the onset of this instability (Balbus 
& Hawley 1991) for a simple case where an axial magnetic 
field Bqz threads a thin Keplerian disc which rotates with 
angular velocity D, = {GAh/r^Y^^ . For axisymmetric per- 
turbations of the disc with (5v = [5vr{z,t),Sv^{z,t),Q\ and 
5B = [5Br{z,t),5B^{z,t),Q\ and for perturbations propor- 
tional to eyqi{}kzZ — iut), one finds the dispersion relation 

1 -1 1/2 

-4+^{k.VA^f , (3) 



where 11.4 = Bo/ -/47rp is the Alfven velocity and = [4fi^ + 
2rQ,dn/dr]^^^ is the radial epicyclic frequency of the disc. In 
order for the perturbation to fit within the vertical extent of 
the disc one needs kzh > 1, where h = Cs/f2 is the half- 
thickness of the disc and Cs is the midplane isothermal sound 
speed in the disc. 

The instability will develop if o;^ < which happens if 
{kvAY < — 2rr2df2/dr. For a Keplerian disc this corresponds 
to (kvAY < 3f2^. Therefore, the above-mentioned condition 
that kzh > I implies that the instability occurs only for va < 
Cs, or the condition for plasma parameter becomes: 



2ci 



Bll 



> 1 



(4) 



where matter pressure p = pc^ . Note that /3 is based on the 
initial vertical magnetic field. Bo- As a result of the instability 
the magnetic field may grow to values much larger than Bo. 
The maximum value of the growth rate is S5(aj)max = 3f2/4, 
and it occurs for fcnmx = (15/16)^^^r2/i>A. For /? < 1 the per- 
turbation does not fit inside the disc and there is stability. As 
P increases from unity (weaker magnetic field) the maximum 
growth rate stays the same but the wavelength of the per- 
turbation gets shorter (oc /3~^''^). For sufficiently small wave- 
lengths the damping due to numerical viscosity (-^ fnumfc^) 



will be larger than the MRI growth rate. To generate MRI- 
driven accretion in numerical simulations, one should have 
large enough /? in the disc so as to have instability, and at the 
same time high enough grid resolution to avoid the damp- 
ing of small wavelengths. This theoretical analysis helps in 
understanding the initial stages of the MRI instability in our 
simulations. 



2.4 Numerical setup 

2.4.1 Method and grid 

We use three-dimensional C3D) MHD second-order Godunov- 
type code developed by our group (Koldoba et al. 2002 ). It 
has many specific features which are oriented towards effi- 
cient calculation of accretion on to a star with a tilted dipole 
or more complex magnetic fields: (1) the magnetic field B 
is decomposed into the "main" dipole component of the star. 
Bo, and the component Bi induced by currents in the disc 
and the corona ( [Tanaka|[l994j ); (2) the MHD equations are 
written in a reference frame which is rotating with the star; 
and (3) the numerical method uses the "cubed sphere" grid 
which has the advantages of both the cartesian and spheri- 
cal coordinate systems and avoids the singularity on the polar 
axis in the spherical coordinate system. The grid on the sur- 
face of the sphere consists of six sectors with the grid in each 
sector topologically equivalent to the grid on a face of a cube 
(e.g., Ronchi et al._1996 ). In contrast with these authors, we 
use a Godunov-tjqae numerical scheme similar to the one de- 
scribed by Powell et al. (1999) and perform simulations in 
the range of azimuthal angles (0 — 2tv) in three dimensional 
space. 

In addition, to better resolve the MRI-driven turbulence 
inside the disc, we compressed the grid towards the equatorial 
plane (see Fig.[T]l. This helps to increase the grid resolution in 
the vertical direction, which is favorable for the investigation 
of MRI-driven accretion and helps to save computing time. 
A number of angular grid resolutions has been used: — 
Ny = 51, 61, 71, 81, 91 (in each of the 6 blocks of the cube). 
MRI-driven turbulence has been observed in all these cases, 
and the picture of the accretion flow is qualitatively the same. 
To save computing time, we use the grid TV^, = A*',, = 61 
as a base. The number of grid cells in the radial direction is 
Nr = 200 (in the strong field case) and A''^ = 220 (in the 
weak field case), so that the region is 1.6 times larger in the 
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case of a weak field. The grid compression has been done in 
such a way that the number of grid cells across the disc in 
the vertical direction is 80, which is sufficient for resolving 
the MRI-driven turbulence. The 3D MHD code is parallelized 
using MPI. We typically use Np = 240 — 360 processors per 
run. 



2.4.2 Reference Units 

3D MHD equations are written in dimensionless form and re- 
sults can be applicable to a wide variety of stars. We chose 
mass of a star as a reference unit for mass, A/o = A-'h and 
radius of the star as reference unit for length, _Ro = ^* ■ Q 
Reference value for velocity is Keplerian velocity at radius 
Ro- vo = [GMo/ RoY^'^ . The reference time scale is period 
of rotation at Rq: Pq — 2itRo/vo. From dimensionalization 
of equations, we get ratio: pova = Bq, where Bo and po 
are the reference magnetic field and density at Rq. We take 
the reference magnetic field Bo such that the reference den- 
sity is typical for considered tjqses of stars (see Tab.[T]l. We 
then define the reference dipole moment po ~ BoRq, density 
Po = Bo/vo, mass accretion rate Mo = povoRo, energy per 
unit time Eo ~ PoVoRo- Temperature To = po/{TZpo), where 
TZ is the gas constant. 

We fix these reference units (and hence, the initial den- 
sity, pressure, etc. distribution in the disc), but vary the initial 
magnetic field on the surface of the star and in the disc with 
dimensionless parameters B'^ and respectively]^ so that 
the dimensional magnetic fields on the surface of the star 
(at the equator) is: B* — B'^Bo and initial field in the disc 
Bd = B'^Bo. 

The MHD equations are solved usingjlimensionless vari- 
ables r — r/Ro, V = v/vo, t = t/Po, B = B/Bq, and so 
on. In the subsequent sections, we show dimensionless val- 
ues for all quantities and drop the tildes. Our dimensionless 
simulations are applicable to stars of different scales. We list 
the reference values for a few types of stars in Tab. [l] To de- 
rive the real, dimensional values, one needs to multiply the 
dimensionless values obtained from simulations by reference 
values presented in Tab. [T] 

2.4.3 Initial conditions 

A rotating magnetized star of mass i\f* and radius i?* is sur- 
rounded by an accretion disc and a corona. The disc is rela- 
tively cold and dense, while the corona is hot and rarefied. 
We set the temperature in the corona = To and in the disc 
Td = COlTc (we determine these values at the inner edge 
of the disc, which is our reference point) . The density in the 
disc and corona pd — po, Pc = O.Olpd. Here, the subscripts 
'd' and 'c' denote the disc and the corona. The initial density 
distribution is derived from the balance between the gravi- 
tational, centrifugal and pressure gradient forces, as well as 

^ The magnetospheric radius, r™. is another important scale of the 
problem, which is one of the main parameters used in theoretical 
investigations. However, in simulations it varies in time, and hence 
we use radius of the star as a reference unit, and obtain rm from 
simulations. 

^ Note that parameter B'^ in this paper is identical to parameter p, 
used in our previous papers (e.g., inlRomanova et al.|2011al. 
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Table 1. Sample reference values for typical CTTSs, white dwarfs, 
and neutron stars. The dimensional values can be obtained by multi- 
plying the dimensionless values by these reference cvalues. 
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Figure 5. Density distribution in the hot spots at the surface of the 
star in the BL regime i B'^ = 0.001) during the local maximum of 
the accretion rate, T = 142 (left panel) and local minimum T = 146 
(right panel) . 

the condition that the corona rotates with the angular veloc- 
ity of the disc and the gas is initially barotropic. In addition, 
the pressure is balanced at the boundary between the disc and 
corona. The resulting initial distributions of density in the disc 
and corona is almost homogeneous (within 10%). The density 
only slightly increases outward in both, disc, and corona. Dis- 
tribution of the temperature (and the sound speed Cs) is also 
almost homogeneous. In dimensionless units, initial values at 
the reference point are: pd = 1, pc = 0.01, Td = 0.01, Tc — 1. 

The star rotates slowly with a fixed angular velocity SI* 
such that the corotation radius r-cor ~ lOi?*, so that fcor /rm ~ 




110 120 130»,„„140 150 160 180 ,. 200 220 



Figure 6. Accretion rate at the surface of the star in cases of the BL 
regime. Panels a and b show the cases of B'^ = 0.001 and B'^ = 1 
respectively. 
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3. We take such a slow rotation to be sure that conditions 
are favorable for magnetospheric accretion during the whole 
simulation run, as we discussed in Sec. |2.2| Cases of more 
rapid rotation and the propeller regime will be investigated 
in the future. 

A star has a dipole magnetic field. The strength at the 
equator is determined by the dimensionless parameter B^. 
The dipole moment is tilted at an angle O about the rotational 
axis of the star (which is aligned with the rotational axis of the 
disc) . Simulations show that at tjqsical parameters of the disc, 
described above, the disc is truncated by the magnetosphere 
(and hence, the magnetosphere is dynamically important) if 

> 1. At smaller values of the field, the magnetosphere 
is dynamically unimportant, and matter accretes in the BL 
regime. We take very weak field, Bi = 0.001, for modelling 
the BL regime. We also take the marginal field, B^ — 1, for 
test runs. A stronger field, B'^ = 10, is used for modelling the 
magnetospheric regime of accretion. 

The disc and corona are threaded with a homogeneous 
small seed poloidal field. The strength of the field is deter- 
mined by the dimensionless parameter B'^. We take — 
0.005 in case of the BL regime, and = 0.01 in case of 
the magnetospheric regime. Figure [2] shows the initial distri- 
bution of density and magnetic field in a;z— slices in cases of 
the BL regime, where Bl = 0.001, 9 = 5° (left panel), and 
magnetospheric regime, where B^ = 10, and the tilt is either 
large, 9 = 30° (middle panel), or very small, 9 = 2° (right 
panel). 

Conversion of our dimensionless parameters to dimen- 
sional shows that the density, velocity, magnetic field and 
other parameters match well with those observed in tjqsical 
stars listed in Tab. ^ However, the temperature in both, disc 
and corona, is about 10-20 times higher than that observed 
in real stars. In our initial setup, this is a consequence of 
the low density contrast between the disc and corona, which 
is initially Pc/pd = 10~^. At the higher density contrast, 
Pc/pd = 10~^, temperature is realistic, however these simu- 
lations require higher grid resolution and the simulation time 
is much longer Realistic, much lower values for temperature 
in the disc and corona are used in axisymmetric simulations 
of a— discs (e.g.,|Lo ng et al. 2005 ) and MRI-driven discs (Ro- 
[manova et al.|2011a) . These simulations do not show any no- 
table difference in results between high and low-temperature 
flows. We suggest that in both cases, the inner disc is suffi- 
ciently thick, so that the gravity force is the main force pulling 
matter to the funnel flow (if < rcor). We discuss possible 
influence of the higher temperature in the disc to magneto- 
spheric accretion in Sec. |4.1| 

2.4.4 Boundary conditions 

At the outer boundary of the simulation region, we fixed all 
8 variables. The boundary is located at _Rout ~ 51i?* for the 
case of the BL regime and _Rout ~ 32_R* for the magneto- 
spheric regime. We observed from simulations that the disc is 
sufficiently large to supply matter during the whole simula- 
tion run. 

At the inner boundary (which is a stellar surface) we 
take 'free' boundary conditions, which are corrected to the 
specifics of the considered problem. In general, the number 
of the boundary conditions should correspond to a number 
of waves propagating from the boundary to the simulation 



region. However, the relationship between the radial veloc- 
ity at the inner boundary and velocities of MHD waves is not 
known in advance. In addition, the interaction of the accret- 
ing matter with the stellar surface is a separate, complex prob- 
lem, which requires detailed calculations of processes at much 
smaller scale Q Besides, matter can accrete, or, it can flow 
away from the star. We do not consider the complex processes 
near the stellar surface, but instead use simplified boundary 
conditions, which allow the absorption of matter which falls 
to the star's surface. 

We consider two types of boundary conditions for dif- 
ferent regimes of accretion. Type A: If the magnetic field of 
the star is dynamically unimportant and accretion is super- 
Alfvenic, then we use 'free' boundary condition for all vari- 
ables: dU/dr — where U = s (entropy), v (total veloc- 
ity vector), Bt (tangential to the stellar surface component 
of the field), r^Br (normal to the surface component of the 
field multiplied by r^). We artificially decrease the pressure 
p at the inner boundary so that to stimulate accretion, and 
to remove excessive matter if it accumulates at the bound- 
ary: po = kpi, k = 0.9, where the subscripts '0' and '1' 
mark the boundary grid and the last calculated cell respec- 
tively. The density in the boundary cell is recalculated as 
Po ~ {Po/so)^^''- These conditions are used in the regime of 
the boundary layer accretion. 

Type B: In case when the magnetic field is dynamically 
important, then the accreting flow is sub-Alfvenic (at the sur- 
face of the star) and two waves propagate from the bound- 
ary to the simulation region: Alfven wave and fast magne- 
tosonic wave. Hence, we provide two additional boundary 
conditions: in the coordinate system rotating with the star, 
the velocity vector is parallel to the magnetic field vector, 
which we split to two components: normal (or, radial) compo- 
nent, Br, and the component tangential to the surface of the 
star, Bt. This condition means that the strong magnetic field 
acts as 'rails' attached to the non-deformed surface of the star 
along which matter flows to the star. We use the following 
algorithm to calculate these conditions: 1) we calculate the 
total vector of velocity in the boundary grid as vq = vi; 2) 
we subtract rotation of the star: Vg = vq — (fi* x ro), where 
ro is the radius-vector of the center of the boundary grid, 3) 
we project vector vq to the direction of the magnetic field: 
Vq = (BoVo)/i3o; 4) we add to the velocity a part con- 
nected with rotation of the star: vq = V(, + (fi* x ro). We use 
these boundary conditions for modelling the magnetospheric 
regime. 

2.5 Number of MRI wavelengths per thickness of the 
disc 

Before starting simulations, we estimate the expected number 
of MRI wavelengths within the thickness of the disc, using the 
analysis outlined in Sec. |2.3| and the dimensionless parame- 
ters of the model. The wavelength of the most unstable, MRI- 
driven mode is Amri ~ 2tvva,z/^k. In the approximation of 
the thin disc, the full thickness of the disc is 2h ^ 2cs/^k, 
where Cs = \/p/p is the locally isothermal sound speed in the 

^ For example, interaction of the supersonic funnel stream with stel- 
lar surface leads to the formation of the radiative shock wave, which 
oscillates (e.g., |Koldoba et al.|2008| . 
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Figure 7. A three-dimensional view of accretion on to a star with a large magnetic field, B'^ = 10, at large © = 30° (left-hand panel) and 
small = 2° tilts. The green background shows one of density levels, the blue color shows the equatorial density slice, and the lines are sample 
magnetic field lines. 




disc. The number of wavelengths per thickness of the disc is 
A^'mri = 2/i/Amri ~ Cs/tvva,z- In our model (in dimension- 
less units) : density in the disc is pd ~ 1, the sound speed is 
Cs ~ 0.1, the Alfven velocity (based on B'^ = 0.01 field) is 
VA,. = B'J^/W^ ^ 2.8 X 10^^(By0.01)(pd/1.0)~^''^ Sub- 
stituting these values to the initial formula for A^mri and to 
eq.[5] we obtain: 

These formulae are approximate but useful for understand- 
ing the start-up conditions for development of the instability. 
Below, we show results of our simulations in two different 
regimes: the BL regime, and regime of the magnetospheric 
accretion. 



3 BOUNDARY LAYER REGIME 

First, we investigate the case where the magnetic field of the 
star is dynamically unimportant and matter accretes to the 
star in the BL regime. We use this case to investigate the MRI- 
driven turbulence in the disc in our numerical setup. This case 
is also interesting from the astrophysical point of view, be- 
cause many stars are expected to accrete in the BL regime. 
Below, we discuss both aspects of the problem. 

3.1 MRI-driven accretion 

We investigated accretion at several values of the seed 
poloidal field in the disc from B'^ = 0.002 up to B'^ = 0.01, 
and found that in all cases the MRI-driven turbulence is de- 
veloped, matter accretes inward, and angular momentum is 
transported outward by the magnetic shear stress (see details 
is Sec. [5|. In experiments with the smallest in the set mag- 
netic field, B'd — 0.002 1^ the initial number of wavelengths 

* We used this field as a base in axisymmetric simulations ( |Ro-| 
|manova et al.|2011a^ . 
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Figure 9. Top panel: The density distribution in the inner disc and in funnel streams in the xz— plane and selected field lines. Bottom panel: Linear 
distributions of different variables in x-direction in the equatorial plane: the density,p, two components of magnetic field, Bz and B^, and 
The dashed vertical lines correspond to /3i 
star. The position of the corotation radius, r 



1 and show the position of the magnetospheric radius to the left, rmi and to the right, 
3r, is also shown. 



, of the 



per thickness of the disc is relatively large, A'^mri = 55 (at 
the inner disc), but the accretion rate, which is proportional 
to the magnetic stress in the disc, is low?: matter of the inner 
disc reaches the star in T ^ 200 rotations (at r — 1). These 
simulations last up to T ^ 800. In the case of the largest 
field, B'^ = 0.01 (A^mri = 11), the accretion rate is high, 
however only 1-2 turbulent cells per thickness of the disc are 
observed in developed turbulence. We take the intermedi- 
ate field, B'^ = 0.005, as a base for simulations, because the 
accretion rate is relatively high, and at the same time, a few 
turbulent cells per thickness of the disc are observed. 

We also varied the magnetic field of the star from very 
weak, Bi = 0.001 (and dynamically insignificant during the 
whole simulation run), to somewhat stronger, iJ^ = 1, where 
a tiny magnetosphere was observed initially, but it disap- 
peared later, when the accretion rate become larger. 

The initial density and temperature distribution in the 
disc is almost homogeneous (within 10%), and the magnetic 
field threading the disc and corona is a constant (see Fig. |2] 
left panel, and Sec. |2.4.3[ l. As a result, the initial distribu- 
tion of /?— parameter in the disc is almost homogeneous, and 
for the disc field of B'^ = 0.005, it is /3 9, 000. We ob- 
served, that the initial number of wavelengths per thickness 
of the disc is A'^mri ~ 10 — 20, as predicted by the theory 

^ Note, that the initial value of A^mri is based on the initial seed 
poloidal field in the disc. Subsequently, the azimuthal component in- 
creases, and the number of turbulent cells per thickness of the disc 
drops. 



(see Sec. |2.5| l. Subsequently, the perturbations grow due to 
the magneto-rotational instability, and turbulence develops in 
the inner parts of the disc (at T ^ 10 — 20 periods of rota- 
tion at r = 1). Later, turbulence develops at larger radii. Sub- 
sequently, the magnetic energy in the disc increases, mainly 
due to the growth of the azimuthal component of the field, 
which leads to strong decrease of /?. As a result, the number 
of turbulent cells per thickness of the disc decreases up to 2-3. 

Figure [5] shows slices of the density distribution at T = 
150. The left-hand panels show that the flow is turbulent, and 
the poloidal field lines are tangled and have a number of re- 
versals which track different turbulent cells. The right-hand 
panel of Fig. [s] shows the equatorial slice of the density dis- 
tribution. One can see that the turbulent cells are strongly 
elongated in the azimuthal direction, and they also look like 
parts of spiral waves. 

3.2 Azimuthally-wrapped magnetic field in the inner 
disc 

In the inner parts of the disc, the azimuthal magnetic field 
increases rapidly due to strong shear in the flow. In addition, 
the poloidal magnetic field increases towards the inner disc 
due to the convergence of the flow and conservation of the 
poloidal magnetic flux. As a result, the azimuthal component 
of the field reaches the value of ^ 1 which is about 200 
times larger than the initial seed poloidal field of B'^ = 0.005. 
The plasma parameter f3 drops to /3 « 1 — 3, that is, the mag- 
netic and thermal pressure are close to the equipartition state. 
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Fig. |4| (right panel) shows the distribution of the mag- 
netic pressure scaled to matter pressure, /3^^ = B^/Sivp at 
t = 150. The thick black line shows exact equipartition, where 
/3 = 1. One can see elongated islands, where the magnetic 
pressure dominates (yellow and red color). The left-hand 
panel shows the 2:2:— slice, where the turbulent cells are seen 
in their cross-sections. These simulations are in accord with 
earlier simulations by Armitage ( 2002) and by Steinacke r &[ 
|Papaloizou ( 2002) , who also observed strong amplification of 
the field in the inner disc. 

|Pringle| l[l989^ argued that magnetic flux built up in the 
disc should escape to the corona due to the Parker instability, 
and that this field can be responsible for magnetically-driven 
outflows from young stars (see also ,Shu et al.]|1988p . For- 
mation of the magnetically-dominated corona has been ob- 
served in 3D simulations (e.g., 'Miller & Stone'2000t |Ha"wley| 
et al.||2 002]). Ver y long-lasting axisymmetric simulations by 
Romano va et al.| ([2011a) performed at 10 times lower den- 
sity in the corona (compared with present simulations) show 
the formation of a very large magnetic corona which slowly 
moves away from the disc to large distances, and which is 
driven by the magnetic force (see fig. 20 and 21 of |Romanova| 
|et al.|[2011a[ ). In the present simulations, we have not seen 
such a corona, because the matter density in the corona is 
too high, and simulations are not as long. We expect that the 
large-scale expanding magnetic corona will be seen in the fu- 
ture 3D simulations. Formation of faster outflows from the 
disc-star, boundary, as discussed by |Pringie| ( |1989[ l, also can 
no be excluded. 

We also observed an interesting phenomenon, that the 
magnetic flux of the star increases in time, and a new mag- 
netosphere forms around the star. The initially weak mag- 
netic field of = 0.001 increased up to B* ^ 1 (at the 
pole, above the star) because the accreting matter brings the 
poloidal field towards the star, and this field is accumulated at 
the stellar surface. Fig. [4] (left panel) shows that the new mag- 
netosphere has a monopole-type shape. The magnetosphere 
does not truncate the disc[^ and the disc accretes in the BL 
regime. We suggest that in longer simulation runs, the mag- 
netic field can increase up to even larger values. On the other 
hand, the disc can bring matter with the field of one or an- 
other polarity, and the induced magnetic field of the star will 
flip its polarity. 

Of, course, details of the disc-star interaction depend on 
the boundary conditions at the star. Here, we show results 
for one type of boundary conditions (Type A, see Sec. |2.4.4l l, 
where both, matter and magnetic field can accrete on to the 
surface of the star However, conditions at the boundary can 
be different. For example, if a star is a perfect conductor, then 
the magnetic flux of the disc will not penetrate through the 
star, but instead it will accumulate at the disc-star boundary. 
Such a condition will be more favorable for formation of the 
magnetic torus, and/or for outflows. 



Note that the field BJ^ si 1 corresponds to the border between 
dynamically active and passive magnetospheres. 
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Figure 10. Distribution of the kinetic energy flux in the hot spots at 
the surface of the star as seen from the top (—2 -direction, left-hand 
panel) and the bottom (+z-direction, right-hand panel). 



3.3 Equatorial belt spots and variability associated with 
accretion of turbulent cells 

In the BL regime, the disc matter accretes to the star in the 
region of the equator, and it forms the belt-shaped hot spot. 
About a half of the gravitational energy of the accreting mat- 
ter is expected to be released in the equatorial belt spot. The 
details of the physics of the disc-star interaction including the 
processes of the heating and cooling of matter in the BL are 
not established. The interaction may include complex pro- 
cesses such as the Kelvin-Helmholtz instability and the for- 
mation of a turbulent layer at the surface of the star,. One can 
expect, however, that most of energy will be released in re- 
gions of enhanced density. This is why we show the density 
distribution at the star's surface. 

We observed in simulations that matter accretes onto the 
star in azimuthally-stretched turbulent cells. Near the star, 
the turbulent cell becomes stretched in the azimuthal direc- 
tion, and when it accretes on to the star, it forms the higher- 
density spot which is also spread in the meridional direction. 
The spots are less dense and more narrow during periods of 
lower accretion rate. 

Fig. [5] (left panel) shows the density spot during accre- 
tion of one of the turbulent cells. The spot is not symmetric, 
because the turbulent cell is not azimuthally-symmetric. The 
spot is quite wide in the meridional direction. The right panel 
shows the spot during period of the minimum of the accre- 
tion rate, when one of cells already accreted, while another 
one did not approach the star yet. This spot is more narrow 
in the meridional direction. 

Accretion of turbulent cells leads to variability of the ac- 
cretion rate at the star Fig. [6] shows the accretion rate for 
our reference case with initial stellar field of = 0.001 (left 
panel), and in the case of larger initial field of the star, — 1 
(right panel) . For these cases, we used boundary conditions of 
Type A and Type B respectively (see Sec. |2.4.4| l. We observed 
that the peaks in the accretion rate (see arrows in the plot) 
correspond to the accretion of individual turbulent cells. The 
accretion rate is typically about 20% higher than the average 
during accretion of a turbulent cell, and the maxima at 50% 
correspond to the accretion of the largest turbulent cells. The 
time interval between peaks in these two cases is ^ 5 and 
At Si 10, respectively. The accretion rate curves are quite sim- 
ilar in spite of the different initial fields of the star and some- 
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what different boundary conditions. Note that we have the 
same initial conditions in discs, and hence the MRI-induced 
turbulence in the disc has the same properties. This may de- 
termine the similarity in the time-scale between peaks. Some 
difference is connected with different conditions at the star. 

In our simulations, we observe quite large turbulent cells, 
and hence we observe and 'count' them during their accretion 
on to the star. Such turbulence dominates in discs which are 
close to equipartition with /9 ~ 1. We suggest, that equiparti- 
tion is a probable state in the inner regions of accretion discs 
where the shear is high. In the opposite case, /3 >> 1, mul- 
tiple turbulent cells are present in the disc. The simultaneous 
accretion of multiple cells will lead to much smoother accre- 
tion, and results will be closer to those observed in cases of 
a— discs (e.g., [Romanova et al.|2004| . Hence, the variability 
pattern carries information about the property of the turbu- 
lence in the disc. 



4 REGIME OF MAGNETOSPHERIC ACCRETION 

Next, we discuss accretion on to stars with dynamically impor- 
tant magnetic fields, where the disc is truncated by the mag- 
netosphere of the star and matter accretes in funnel streams. 
We chose the parameter — 10 for these runs. We inves- 
tigate cases of large (9 = 30°) and small (0 = 2°) tilts of 
the dipole magnetic moment relative to the rotational axis. 
To save computing time, we increased the seed magnetic field 
in the disc up to = 0.01. This helped to increase the accre- 
tion rate (which is proportional to the magnetic stress in the 
disc), and cause the disc to reach the surface of the star more 
rapidly at about T = 4(p^ 

We observed that initially, about 10 MRI wavelengths fit 
within the vertical extent of the disc, and the inner parts of the 
disc started accreting towards the star. Later, the azimuthal 
field increased, and the number of turbulent cells observed 
later in simulations dropped to 1-2 per thickness of the disc. 

Figure[7]shows a 3D view of the MRI-driven discs in cases 
of low and high tilts of the dipole field. The background shows 
one of density levels. One can see that in both cases, the disc 
is strongly inhomogeneous: matter tends to be concentrated 
in azimuthally-elongated spiral segments. The magnetic field 
lines are stretched in the azimuthal direction. Below, we con- 
sider these two cases in greater detail. 

4.1 Accretion on to a star with a large tilt, 6 = 30° 

In the case of a large tilt of the dipole field, matter flows to 
the star in funnel streams. Figure [t] (left panel) shows such an 
accretion flow, where matter of the disc accretes to the star 
in funnel streams. These funnel streams form when a large 
turbulent cell moves towards the magnetosphere, where it is 
stretched, slowed down by the magnetosphere, and pulled to- 
wards the star by the gravity force. 

Figure [s] (left two panels) shows x — z and y — z slices 
of density distribution and the poloidal magnetic field. One 
can see that the disc is truncated and matter flows to the star 
in funnel streams. The height of the magnetosphere is only 

Simulations with a field of B'^ = 0.005 produce similar results, 
but on a longer time scale. 



slightly larger than the height of the turbulent disc. The re- 
gions of enhanced density in the disc (yellow and red colors) 
track the turbulent cells. We observed accretion of several tur- 
bulent cells during the simulations. 

The field lines threading the disc inflate into the corona. 
The accretion disc compresses the external field lines and 
pushes them to reconnect. After reconnection, the magnetic 
islands form and propagate into the corona. This process of 
compression and reconnection is similar to that observed in 
axisymmetric simulations ( [Romanova et al.||2011a ). Recon- 
nection acts as an efficient means of diffusivity, and it helps 
the disc to penetrate through the external layers of the mag- 
netosphere of the star towards the closed magnetosphere and 
to the star through the funnel flow. 

The right-hand panel of Fig. [8] shows that the magnetic 
field has a significant azimuthal component, and the turbu- 
lent cells are elongated in the azimuthal direction. The tilted 
rotating dipole excites two spiral waves in the inner disc. 
These waves are similar to those observed in 3D simulations 
of a— discs with a tilted dipole magnetic field ( [Romanova et| 
|al.|2011bl ). 

Figure [9] (top panel) shows the disc-magnetosphere 
boundary in greater detail. One can see that the disc is tur- 
bulent. It is disrupted by the stellar magnetosphere and the 
matter flows towards the star forming a funnel flow. One can 
see that the funnel streams are not symmetric: the accretion 
rate is higher from the left side than from the right side. This 
is in contrast with earlier 3D simulations of accretion from 
a— discs, where the funnel streams are symmetric (see, e.g.. 
Fig. 4 from [Romanova et al.[[20(34| . The asymmetry is con- 
nected with the fact that matter approaches the magneto- 
sphere of the star and then accretes in large-scale turbulent 
cells, which are usually closer to one of poles than to another. 
Fig. [t] (left panel) shows a tjqsical process of accretion of the 
turbulent cell. 

The bottom panel of Fig. [9] shows the distribution of dif- 
ferent variables in the equatorial plane. One can see that the 
density p is highly variable in the inner disc, and it is very 
small inside the magnetically-dominated magnetosphere. The 
azimuthal component of the field, B^, is a few times larger 
than iJz -component. Inside the magnetosphere, the toroidal 
field is much smaller than the poloidal field. 

The inner radius of the disc coincides with the point 
where our modified plasma parameter parameter /3i (eq. 2) is 
unity. The vertical dashed lines in Fig. [8] show the radii where 
Pi — 1. One can see that the magnetospheric radius to the left 
of the star, rmi ~ 2.4 is smaller than to the right, rm2 ~ 3. 
The asymmetric nature of the magnetospheric surface is also 
seen in Fig. [8] (right panel). Top panel of the Figure also shows 
that the funnel streams have a finite width and some matter 
start flowing to the funnel from larger radii, r ^ 3 — 4. 

Matter moves to the funnel flow due to the gravitational 
force, which dominates over the centrifugal force (in our case 
of a slowly rotating star). The lifting of matter to \z\ > to 
the funnel starts at the radius where the azimuthal motion 
of the inner disc matter slows down due to interaction with 
the magnetosphere. This occurs at r ^ 4 — 5. The magneto- 
sphere is small, while the disc is has a finite thickness due to 
the turbulent nature of the flow. Hence the magnetosphere is 
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not a big obstacle for the flow of matter towards the staip^ 
The matter pressure gradiem force also contributes to the lift- 
ing, though the role of this force is not as significant as in 
the case of the aligned dipole (e.g., |Romanova eFal.|20 02 and 
|Campbell|2010, ) Note, that axisymmetric simulations of MRI- 
driven accretion were performed at much lower, realistic tem- 
perature in the disc and corona, and did not show principal 
difference between higher and lower temperature cases {Ro-| 
[manova et al.|2011a[ ). Moreover, in case of a tilted dipole, the 
magnetic field lines are inclined in the equatorial plane, and 
hence the star's gravitational force can pull matter directly 
from the equatorial plane of the disc. 

Note, that inside the disc, the total matter stress is much 
larger than the magnetic stress, /3i >> 1, and the disc acts as 
matter-dominated (see also Sec.[5|. The MRI-driven instabil- 
ity, initiated by the initial seed field in the disc, provides in- 
ward transport of matter, which can be roughly described as 
an effective 'viscosity'. However, it does not provide diffusivity 
at the disc-magnetosphere boundary, he magnetic field of the 
star does not thread the incoming disc and does not influence 
to the MRI processes inside the disc. As a result the matter- 
dominated disc compresses the magnetosphere and pushes it 
to reconnect as shown in Fig. [s] and Fig.|9] It also penetrates 
toward the closed magnetosphere due to instabilities. 

The funnel stream hits the surface of the star with a high 
velocity and the kinetic energy of the flow is released. Figure 
[To] shows the distribution of the kinetic energy flux in the hot 
spots on the surface of the star due to the impact of the funnel 
streams. One can see that the spots forming at the south and 
north hemispheres of the star are different. The difference 
between spots reflects the difference between funnel streams 
approaching the south and north hemispheres. Note, that hot 
spots obtained in 3D simulations of a— discs, are almost iden- 
tical ( [Romanova et al.|2004[ l. Inequality of spots in the case of 
turbulent discs may have important consequences for analysis 
of variability from accreting magnetized stars. 



4.2 Accretion on to a star with small tilt, 6 = 2° and 
interchange instability 

Here, we discuss results for the case of very small tilt, = 2°. 
Figure [TT] (left panels) shows that the magnetosphere inflates 
in an almost symmetrical fashion. The magnetic flux inflates 
sideways, as we observed in axisjmmetric simulations (Ro^ 
[manova et al.|[2011a[ ). This is because the time scale of disc 
accretion is smaller than the inflation time scale. The mag- 
netic field in the inflated and stretched corona reconnects, 
while magnetic islands form and propagate outward. 

The right-hand panel of Fig.[Tl]shows an equatorial slice. 
One can see that the magnetosphere with such a small tilt 
does not excite strong spiral waves. However, we observed 
another interesting phenomenon: matter of the disc located 
at large distances from the star penetrates through the mag- 
netic field of the star in thin dense filaments. This is a sign 
of interchange (or magnetic Rayleigh-Taylor) instability. Fig. 
[t] (right panel) also shows these filaments. The filaments are 
tall (in the vertical direction) and narrow (in the radial di- 
rection) . This shape is favorable for penetration between the 
field lines. This situation is analogous to penetration of disc 
matter through a vertical magnetic field threading the disc. 
This was investigated earlier with analjrtic theory and MHD 
simulations (Kaisig et al.|1992[ [Spruit et al. 1995j|Lubow &| 
[Spruit[[T995[ l. Here, we observe this phenomenon in simu- 
lationsp'l The interchange instability can be an important 
mechanism of matter penetration through the magnetic field 
in accretion discs in different situations. It can give an en- 
hanced transport of magnetic flux. 

Matter of the disc which reaches the disc-magnetosphere 
boundary, can further penetrate through the magnetically- 
dominated magnetosphere due to same instability (e.g., Arons[ 



& Lea 1980; Spruit & Taam 1990} |Li & Na rayan 20 041) or 
due to the Kelvin-Helmholtz instability (e.g., Lovelace et al.[ 
[2010p . Our earlier global simulations of a— type discs have 



1^ Note, that the situation is different in cases of much large magne- 
tospheres of X-ray pulsars and strongly magnetized white dwarfs in 
cataclysmic variables, where the hight of the magnetosphere may be 
orders of magnitude larger than the thickness of the disc. 



^■^ Such an instability has been also observed in global 3D MHD 
simulations by I gumenshchev[ l|2008), where the magnetic flux was 
accumulated in the inner parts of the disc around the black hole (see 
also |Narayan et al.|2003p . 
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shown that in magnetized stars matter flows towards the stel- 
lar surface either through two ordered funnel streams, above 
the magnetosphere, where two hot spots form, and the light 
curve is almost sinusoidal. Or, it can penetrate the magneto- 
sphere due to interchange instability in the equatorial plane 
in narrow and tall 'tongues'. In this case, a number of tempo- 
rary hot spots form in random places, and the light-curve is 
chaotic C Romanova et a l.'2008VK ulkarni & Romanov a|2008p . 
This has an important application for understanding of peri- 
odic and quasi-periodic variability in differen stars with small 
magnetospheres (e.g., |Kulkarni & Romanova||2009l |Bachetti| 
let al.|2010D . 

There are several factors which determine the bound- 
ary between the stable and unstable regimes of accretion. 
The unstable regime is favorable, if the magnetic axis has a 
small tilt about the rotational axis, if a star rotates slowly 
(hence, the gravitational plus centrifugal potential is nega- 
tive), and also if sufficient amount of matter is accumulated 
at the disc-magnetosphere boundary. More, precisely, the the- 
oretical condition for instability which agrees well with re- 
sults of simulations requires that the gradient of the surface 
density S per unit of the vertical magnetic field should be 
large and positive, d(T./B^)/ dr > (e.g., [Spruit et al.|199"5 



|Kulkarni & Romanova||2008l ). This term is required to over 
come the stabilizing effect of the velocity shear (see, e.g., eq. 
59 of |Spruitetal.|1995j ). 

In the present simulations, we see that matter is lifted 
above the magnetosphere and accretes in funnel streams (see 
Fig. left panels). However, this is the case, which is fa- 
vorable for the accretion though instabilities: the star rotates 
slowly, and the tilt O = 2° is very small. It is possible however 
that a high disc temperature (see Sec. |2.4.3p drives matter to 
the funnel flow and hence opposes to the onset of the instabil- 
ity. However, this is not the case, because in the a— disc mod- 
els (e.g.,l Romanova et al.|2008 )l, where instability is observed, 
we used the same parameters for the density and temperature 
distribution in the disc and corona. We suggest the instability 
is not seen in the current simulations, because the runs are 
not long enough. There is not enough time to build up the ra- 
dial gradient of (E/Bz) needed to trigger the instability. Note 
also that conditions for the onset of the interchange instability 
may be different in case of a turbulent disc. 



5 ANALYSIS OF STRESSES IN THE DISC 

In both cases considered above the turbulent accretion is ex- 
cited by the magneto-rotational instability. Here, we analyze 
stresses in the disc. We separate the disc from the low-density 
corona using a specific density level, pdiac ~ 0.3, which is 
tjqjical for the boundary between the disc and corona. We in- 
tegrate the matter (subscript 'm') and magnetic (subscript 'f ') 
stresses across the disc (in the z— direction) and obtain 



(Tf) = - 



2h 



dz 



Br B^ 

4-K 



Here 



are averaged azimuthal velocity and matter flux, and h — 
h{r) is the half-thickness of the disc. Note that r{Tm + Tf) is 
the radial angular momentum flux in the disc. 
The matter and magnetic pressures are 



2h 



dzP, {Pf} = 



2h 



dz 



'Sir' 



The standard a— parameters associated with matter and mag- 
netic stresses are: 



2 (T^) 

3 (Pm) ' 



a/ = 



2 (Tf) 

3 (Pn) ' 



Figure [12] (left-hand four panels) shows the equatorial 
distribution of different averaged values in the disc for the 
case of the boundary layer accretion, where the magnetic 
field is relatively weak (Bi — 0.001). One can see that the 
magnetic stress, (T/), is larger than matter stress, (Tm), and 
hence angular momentum is transported by the magnetic 
stress. Note that the matter stress also results from magnetic 
turbulence. It is interesting that the magnetic stress forms 
a clear one-armed spiral wave. The same wave is seen as a 
low-density pattern in Fig. |3] Slices for a— parameters (a/ 
and Om) repeat the pattern for stresses. The four right-hand 
panels of Fig. [12] show the distribution of stresses along the 
X and y axes in the equatorial plane. One can see that the 
magnetic stress is larger on average than the matter stress in 
both the X and y directions. Its distribution is very inhomoge- 
neous and it reflects the distribution of matter and magnetic 
stresses. The a— parameters are also strongly inhomogeneous, 
and they track the regions of enhanced stresses. The maxi- 
mum values are high: q/ ^ 0.2 — 0.35 and Om ~ 0.05 — 0.18. 

Figure [13] (left-hand four panels) shows the equatorial 
distribution of different averaged values in the disc in the case 
of the magnetospheric accretion where the magnetic field is 
relatively strong and the tilt of the magnetic moment is high: 
-B^ = 10, = 30°. Here, we see that the equatorial distribu- 
tion of stresses and a— parameters has a similar pattern. The 
magnetic stress does not repeat the clear two-armed pattern 
observed in the density distribution (see Fig. [8| . 

The matter pressure is much higher than the magnetic 
pressure (because the temperature in the disc is relatively 
high due to the initial conditions). However, the magnetic 
stress is larger than matter stress at most of radii. 

The bottom left-hand panels show that the magnetic and 
matter a— parameters (a/ and a™., respectively) repeat the 
pattern of the magnetic and matter stresses. The bottom right- 
hand panels show that the a/ parameter is largest in spiral 
waves, (q/ ^ 0.1 — 0.3), and it is much smaller between 
them and at larger distances. The matter a— parameter, a™, 
is much smaller than a f . 

Note that the simulations in the case of a weak stellar 
magnetic field, B'^, — 0.001, run longer than in the case of 
a strong field. This may explain the formation of clear spi- 
ral features in spite of the small tilt of the magnetic moment, 
O = 5°. In the future, longer simulations runs may show fur- 
ther growth of the magnetic energy in the disc and more de- 
veloped spiral structure. 



6 CONCLUSIONS AND DISCUSSION 

We performed global three-dimensional ideal MHD simula- 
tions of MRI-driven accretion on to magnetized stars with the 
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Figure 12. The integrated stresses, pressure and q— parameters associated with matter ((T„i), (Pm), ctm) and the magnetic field ((T)), (Pf), 
o in the model with Bl = 0.001 and © = 5° . Left two panels: the equatorial distribution. Right two panels: distribution along the x and y axes. 
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Figure 13. Same as in Fig.n2^but for the case with B'^ = 10, = 30° 



dipole magnetic fields tilted at an angle O relative to the stel- 
lar rotational axis (which coincides with the rotational axis 
of the disc). We use the version of our cubed sphere code 
fKoldoba et a l.|2002| l where the grid is compressed near the 
equatorial plane. Simulations were performed in dimension- 
less form and can be applied to different types of stars with 
relatively small magnetospheres, such as CTTSs, accreting 
brown dwarfs, millisecond pulsars, and dwarf novae. 

Multiple simulation runs were performed for different 
parameters of the star and the disc. However, for clarity we 



show results for two main regimes of accretion. In one of 
them, the magnetic field of the star is djmamically unimpor- 
tant: it does not stop the disc, and matter accretes to the sur- 
face of the star in the BL regime. In another regime, the mag- 
netic field is dynamically important, it truncates the disc, and 
magnetospheric accretion is observed. 

Common features observed in both regimes of accretion 
are the following: 

(i) We observed that the MRI-driven turbulence develops 
in the disc and it is similar to that observed in cases of accre- 
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tion on to non-magnetized objects (black holes) (e.g., |Hawley| 
|et al.|200T] l . The magnetic stress in the disc is larger than the 
turbulent part of matter stress, and hence, the angular mo- 
mentum is transported outward by the magnetic stress. 

(ii) The disc is strongly inhomogeneous: it consists of mul- 
tiple azimuthally-elongated turbulent cells. The azimuthal 
field is about 10 times larger than the poloidal field. 

(iii) The spiral structure is observed in the density and 
stress distributions in the disc. It is more prominent in the 
case of strongly tilted, dynamically important field (B* = 10, 
6 = 30°), and in case of the dynamically unimportant field 
(B^ — 0.001, 6 = 5°) where simulations are longer. 

Below, we show main results for our two regimes of accre- 
tion: 

Boundary Layer regime: 

(iv) The magnetic field of the inner disc is strongly ampli- 
fied due to the shear. The magnetic pressure becomes compa- 
rable with matter pressure, p ^ 1 — 3, and hence the inner 
disc is almost in the equipartition. This result is in accord with 
earlier simulations by Armitage (2002) and Steinacker & Pa-' 
Ipaloizou] l [200 2 ), who also observed strong amplification of 
the field in the inner disc. 

(v) Matter of the disc accretes to the star in individual tur- 
bulent cells. Accretion of the turbulent cell leads to formation 
of the density 'spot' on the star which has a shape of the equa- 
torial belt. The belt is inhomogeneous and it is wide in the 
meridional direction. Spots have lower density and are more 
narrow during periods of lower accretion rate. 

(vi) The matter flux at the surface of the star varies, and 
the maxima in the variability curve correspond to accretion of 
individual turbulent cells. The process is quasi-periodic with 
a period of a few Keplerian rotation periods at the radius of 
the star. 

(vii) The accreting matter brings the poloidal magnetic 
flux to the star, and a weakly magnetized star acquires a new, 
'induced' magnetosphere. The polarity of the induced field de- 
pends on the polarity of the disc field, and hence the field of 
the star can flip its polarity. 

In the Magnetospheric regime: 

(viii) The disc is truncated by the magnetosphere of the 
star at a few stellar radii, where the magnetic stress in the 
magnetosphere balances the total matter stress in the disc. 
Closer to the star matter flows to the star's surface through 
the funnel streams. 

(ix) The funnel streams flowing towards the south and 
north hemispheres are different, because the accreting tur- 
bulent cell flows mainly to the closest magnetic pole, and less 
matter flows to another pole. 

(x) The hot spots on the star are also different. They reflect 
the energy distribution in funnel streams. 

(xi) In the case of a small tilt, = 2°, matter of the disc 
penetrates through the external magnetosphere of the star 
due to the interchange instability. We suggest that this insta- 
bility can also play an important role in further penetration of 
the disc matter through the closed, magnetically-dominated 



In the opposite case of multiple turbulent cells in the disc (not 
considered here), the number of cells will accrete simultaneously and 
the funnel streams will be more symmetric. 



magnetosphere. Further work is required to investigate this 
process. 

One of important feature of accretion from turbulent 
discs (compared with a— discs) is that the turbulence breakes 
the symmetry between two funnel streams and hot spots 
in case of magnetospheric accretion, and it brings non- 
axisymmetry to the shape of the equatorial belt in case of 
the boundary layer accretion. This property is important for 
understanding the light-curves from magnetized and non- 
magnetized stars. 

The further enhancement of the magnetic field at the 
disc-star boundary, may possibly lead to the formation of a 
magnetic equatorial torus ( [Paczyii ski|| 1 9 78| l which has been 
used in the phenomenological models of quasi-periodic os- 
cillation in dwarf novae cataclysmic variables (Warner 2003] 
|Preto rius 2006). On the other hand, the build up of magnetic 
field may flow in to the corona due to the magnetic pressure 
force (see, e.g.. Mille r & Stone |2000[ [ Romanova et al.|2011a| | 
and it may be responsible for outflows from young stars 
( [Pringle|1989, see also Shu et al. 1 988 ). Recent axisymmetric 
simulations show formation of powerful, magnetically-driven 
outflows from the disc-magnetosphere boundary, where the 
strong magnetic field of the star is compressed by the disc at 
very high accretion rate ( jLii et al. i 2011 ). These simulations 
can possibly explain the FU Ori winds as discussed by |K6nigI| 
let al.| l [201ip . We have not seen formation of fast outflows in 
current simulations of accretion in the BL regime. Recently, 
we were able to obtain outflows from stars rotating in the 
'propeller' regime (Ustjoigova et al. 2011). However, this is 
the case of the dynamically important field. 

An interesting phenomenon of magnetospheric accre- 
tion through interchange instability observed in modelling of 
Q— discs (e.g., Romanova et al. 2008; Kulka rni & Romanova| 
2008 ), has not been observed in current simulations. We sug- 
gest that our simulations of magnetospheric accretion were 
not long enough to observe this instability, and/or the insta- 
bility criterion can be somewhat different for cases of turbu- 
lent and a— type discs. 
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